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We present a theoretical treatment of the surprisingly large damping observed recently in one- 
dimensional Bose-Einstein atomic condensates in optical lattices. We show that time-dependent 
Hartree-Fock-Bogoliubov (HFB) calculations can describe qualitatively the main features of the 
damping observed over a range of lattice depths. We also derive a formula of the fluctuation- 
dissipation type for the damping, based on a picture in which the coherent motion of the condensate 
atoms is disrupted as they try to flow through the random local potential created by the irregular 
motion of noncondensate atoms. We expect this irregular motion to result from the well-known 
dynamical instability exhibited by the mean-field theory for these systems. When parameters for the 
characteristic strength and correlation times of the fluctuations, obtained from the HFB calculations, 
are substituted in the damping formula, we find very good agreement with the experimentally- 
observed damping, as long as the lattice is shallow enough for the fraction of atoms in the Mott 
insulator phase to be negligible. We also include, for completeness, the results of other calculations 
based on the Gutzwiller ansatz, which appear to work better for the deeper lattices. 



I. INTRODUCTION 

The transport properties of atomic Bose-Einstein con- 
densates have recently been the subject of much inter- 
est. In a pure harmonic trap, the dipole mode of the 
motion — where the cloud of atoms oscillates back and 
forth without altering its shape — is known to be stable. 
On the other hand, if an optical lattice is used to create 
a one-dimensional array of potential wells and barriers, 
one may find, even in a single-particle picture, a damp- 
ing of the oscillations due to the non-quadratic nature of 
the resulting dispersion relation 0, 0, OJ 0, [M |(| . When 
interactions between atoms are included, at the mean- 
field level, one finds dynamical instabilities Q that may 
result in a very large damping 0, . All these effects 
are, however, only expected to be substantial when the 
quasimomentum of the cloud of atoms is sufficiently large 
(typically, of the order oi irh/A, where A/2 is the lattice 
spacing) . 

In recent experiments with 87 Rb atoms 0, [n| , con- 
fined to move in one-dimensional "tubes," a surprisingly 
large damping of the dipole mode was observed, for very 
weak optical lattices and very small cloud displacements. 
We note that no (or very little) damping was observed 
for the same system in the absence of the tight transverse 
confinement [3 El- In the experiments [T3|, the oscilla- 
tion frequency in the harmonic trap, ujq/2-k was about 60 
Hz, whereas the photon recoil energy Er = h 2 /(2m\ 2 ) 
corresponded to a frequency Ejt/h = 3.47 kHz. Under 
these conditions, for a shallow lattice, the maximum dis- 
placement of the condensate in the experiment (7 to 8 
lattice sites) should not result in a momentum larger than 
about 0.1(7r?i/A), which is well within the quadratic part 
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of the lattice dispersion curve. Likewise, the quasimo- 
mentum spread arising from the finite size of the cloud 
itself was also quite small (of the order of 27r?i/13A, since 
the Thomas-Fermi radius of the cloud is about 13A). 

Since these results were first presented (and, in some 
cases, predating them), a number of theoretical treat- 
ments have been put forward that, directly or indirectly, 
address various relevant aspects of the underlying dy- 
namics, from different perspectives. It has been shown, 
for instance 0, ITHl Ej j that the momentum cutoff for 
the dynamical instability may be substantially lowered 
for commensurate lattices, and, probably more relevant 
for the experimental situation, that the boundary be- 
tween regular and irregular motion becomes "smeared 
out" due to quantum fluctuations. Numerical calcula- 
tions based on a truncated Wigner representation [l7j 
have also shown that the fraction of atoms with momenta 
in the unstable region can indeed cause damping of the 
center of mass motion of the whole system. As we shall 
show below, this fraction is, in fact, a non-negligible num- 
ber, for the experimental parameters, even for relatively 
shallow lattices, because of the large depletion caused by 
the very tight transverse confinement. 

In a recent series of papers 0, 0| , several of us have 
characterized the damping mechanisms that may dom- 
inate, for these systems, in different parameter ranges. 
Perhaps the most important conclusion of these papers 
is that the very deep lattices (lattice potential V larger 
than about 5-Er) can be described very well by an ex- 
tended fermionization model, in which most atoms local- 
ize in a Mott-insulator state with unit filling of the lattice, 
and the remaining atoms are free to move above the Mott 
state with a renormalized kinetic energy. Both the atoms 
in the Mott state and the remaining atoms are treated 
as effective non-interacting fermions whose dynamics are 
governed by a combination of the trap potential and ap- 
propriate kinetic energy terms. These references also 
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show, however, that there is a region of values of the ra- 
tio of interaction energy to kinetic energy (referred to as 
the "intermediate region" in Q ) where the single-particle 
models, whether bosonic or fermonic, are inadequate to 
describe the dynamics of the Bose-Hubbard model, which 
is the main underlying theoretical tool for most of the 
studies described above. This intermediate region, in the 
experiments of ^} > covers all the lattices studied with V 
smaller than about 5Er, although there is some concern 
that for the shallowest lattices the tight-binding approx- 
imation leading to the Bose-Hubbard model itself may 
not be entirely accurate. 



The present paper is an attempt to fill in this gap by 
presenting mean-field based calculations for the Bose- 
Hubbard dynamics in the "intermediate region" where 
the fraction of atoms having undergone the transition to 
the Mott insulator state is still negligible, and the system 
is mostly superfluid, yet the interaction energy cannot be 
neglected. Our main calculational tool is time-dependent 
Hartree-Fock-Bogoliubov theory, and we show that this 
approach docs indeed reproduce qualitatively many of 
the features of the damping observed in the experiments, 
although it generally underestimates its magnitude for 
any given lattice depth. With the insights gained from 
these simulations, we develop a model for the damp- 
ing that leads to a formula of the fluctuation-dissipation 
type, which relates the damping of the center of mass 
motion to the random density fluctuations in the noncon- 
densate atoms that result from the dynamical instability. 
We find that this formula leads to good agreement with 
the experimental data when parameters for the charac- 
teristic strength and correlation times of the fluctuations, 
obtained from the HFB calculations, are substituted in 
it. 



The basic ingredients of our model are: (1) a large non- 
condensate fraction that arises as a direct consequence 
of the enhanced effective on-site interaction (due to the 
tight transverse confinement), of which a non-negligible 
part occupies high-momentum states and is therefore af- 
fected by the dynamical instabilities, and (2) the inter- 
action between these noncondensate atoms and the con- 
densate, which is modeled by exploiting a formal anal- 
ogy with an external, random potential. The latter is 
well known to determine localization of the atomic wave- 
function in one-dimensional systems |l9l l20j| . These in- 
gredients (1) and (2) are introduced in Sections II and 
III, respectively. Section IV then presents the results of 
Hartree-Fock-Bogoliubov (HFB) calculations, which we 
use to estimate the parameters appearing in the damping 
formula. Results from an alternative mean-field theory, 
based on the Gutzwillcr ansatz, are presented in Section 
V. Finally, Section VI is devoted to further discussions 
and conclusions. 
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FIG. 1: Black dots: fractional quantum depletion (h/n = 
1 — n c /n) vs. lattice depth V. Open circles: the fraction of 
n with momenta > nTi/\. Solid line: Eq. © with n — 2.2 
(approximate atom density in the experiment at the center of 
the trap) and a — 0.37 (best fit). 



II. HAMILTONIAN AND STATIC 
(GROUND-STATE) RESULTS 

The starting point for our theoretical treatment is 
a Hamiltonian of the "tight binding" or Bose-Hubbard 
form [HI, 



H = -jJ2 a]&i +SlJ2j% + f - 1) (!) 



<J,t> 



In this expression, the sum < i,j > is taken over nearest 
neighbors, aj(cij) are bosonic field operators that anni- 
hilate (create) an atom at the lattice site j, hj = aja,,-, 
and il = mc^Q A 2 / (8 En) characterizes the strength of the 
harmonic trap. The on-site interaction energy is 



U = 



2ah 



2nE R 



1/3 



(2) 



where Cj x ^ v>z are the oscillation frequencies at individ- 
ual lattice sites along the three axes, obtained with a 
harmonic approximation expanding around the minima 
of the potential wells, a x ,y,z = y/Ti/mCbx^y^ are the ef- 
fective harmonic oscillator lengths and a = 5.31 x 1CP 9 
m is the s-wave scattering length. For the experiment 
Cj x /2tt = Cjy/2-K = 35 kHz and Cj z = \J 4E R s/h, where s 
is the longitudinal periodic potential strength in units of 
the recoil energy, s = V/Er. The "hopping" energy J is 
1 /4 of the bandwidth, and obtained by the usual Mathicu 
function treatment (see for more details). The defi- 
nition (2) above implies that the parameters U, J and £1 
in Q arc understood to be in units of En. The summa- 
tion indices in Eq. (JJJ range from — M/2 to M/2, where 
71/ + 1 is the total number of wells (100-200 in our nu- 
merical calculations), and j = at the center of the trap 
We have used a numerical quantum Monte-Carlo [22| 
method to derive the single-particle density matrix for 
the Hamiltonian at very low temperature (0.01J, in 
our calculations), for a total number of atoms N = 80. 



3 



From it wc obtain the quantum depletion shown in Figure 
1. The curve is a fit to the formula 



operator is 



n 



a 



(3) 



which gives the depiction in the homogeneous case (that 
is, in the absence of the harmonic trap, Q — 0). n c is the 
density of condensate atoms, n the total density, and a 
a parameter that can be calculated from the Bogoliubov 
spectrum of excitations Since, in the trap poten- 

tial, the spectrum is modified and n is not uniform, a in 
the figure has been treated as an adjustable parameter. 
We find that, even for very shallow lattices, some 20 to 40 
percent of the atoms are not part of the condensate. The 
figure also shows (open circles) the fraction of these non- 
condensed atoms that have quasimomcnta greater than 
?nt/\ (also calculated from the numerical single-particle 
density matrix). 

In other studies, we have observed that the Mott in- 
sulator begins to form around V = 3-Er in this system, 
as characterized by a small decrease in the density fluc- 
tuations around the center of the trap that first becomes 
visible at this point. Nonetheless, Figure 1 shows that 
the Bogoliubov result for h remains approximately 
valid until around V = 5-Er, which inspires us some 
confidence that the mean-field analysis that follows may 
be at least semi-quantitatively valid even for those very 
highly-depleted systems. 



III. A DAMPING MODEL 

In some recent work, one of us [25| has developed 
a formalism to describe the effect on matter waves of 
coherence-breaking processes such as random "localiz- 
ing" events, momentum kicks, or perturbation by (time- 
dependent) random external potentials. All these pro- 
cesses can be shown to lead to a damping of the center 
of mass motion of the system. 

Given the relatively large fraction of atoms in the dy- 
namical instability region, calculated in the previous sec- 
tion, it seems, therefore, natural to consider their density 
fluctuations as providing a sort of random potential for 
the condensate to flow through, and to expect the damp- 
ing to arise as a consequence of this. The goal of this 
section is to make this picture plausible and quantita- 
tive, first by redcriving the damping due to an external 
random potential, then by showing how the Heisenberg 
equations of motion derived from the Hamiltonian 1]TJ) 
can be cast into a similar form through a standard fac- 
torization ansatz, and, finally, deriving from all of this a 
damping formula. Results from HFB calculations (to be 
discussed in much more detail in the next section) will 
also be introduced to establish the existence of the requi- 
site random density fluctuations in these systems in the 
relevant regime. 

In our tight-binding model, the center of mass position 



-T 

2N ^ 



(4) 



The commutator of x cm with the Hamiltonian Q yields 
a center of mass velocity operator 



i,J\ m-^, 
2Nh 



](a]a j+1 - a] +1 aj), 



(5) 



which is essentially the same as the current operator in 
pij , A further commutation yields the time derivative 
of a]a J+ i: 

ih^a^aj+i = (2j + l)f2a]% + i + C/a](n J+ i - 



ij+i 



J(a)_ 1 a 



3+1 _ a j a j+2 



T-j+l 



hi) (6) 



with an analogous result for the derivative of Oj +1 cij. 
In the "intermediate region" in which we are interested 
here, where the evolution of the system is not adequately 
described by single-particle models (bosonic or pseudo- 
fermionic), we expect the damping to arise from the in- 
teraction term (proportional to U) in |j(JJ). 

Before we get to work on that term, however, con- 
sider what would happen if one were to replace it in (1) 
by a random external potential, proportional to ^ Vjhi. 
Eqs. (4) and (5) would be unchanged, whereas Eq. © 
would become 

ih^-a^a j+ i =(2j + l)f2a]%+i + (V j+1 - V,)a]a i+1 

+ J(a]_ 1 % + i - ajaj+2 + n j+1 - hi) (7) 

Now consider formally taking the ordinary quantum- 
mechanical expectation value of Eq. J7J) and integrating 
it over an interval (t — At, t), to get 



{a}a j+1 (t)) =(a}a j+1 {t - At)) 



t~At 



(V J+1 -V 3 )(t>)(a% +1 (t>))dt' + 



(8) 



where . . . represents terms that do not contain Vj 's. Sub- 
stituting JHJ back into the (expectation value of the) sec- 
ond term on the right-hand side of JJJ), one obtains two 
kinds of terms: some linear in the Vj, and some quadratic 
in Vj. The linear ones involve products of a Vj at the 
time t and field operators at an earlier time, and we may 
assume that they vanish in an ensemble average over dif- 
ferent realizations of the random process Vj, provided it 
has a sufficiently short correlation time. The ensemble 
average of the quadratic terms, on the other hand, yields 



t-At 



m+i - Vj)(t)(V j+1 - Vj)(t')}(ata j+1 (t'))dt' 



--T c ((V J+1 -V,) 2 )(a]a, +1 (t)) 



(9) 
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where r c is the characteristic correlation time for the 
randomly-fluctuating potential, {{Vj+i — Vj) 2 ) is the av- 
erage (squared) strength of the fluctuations, and it has 
been assumed that (otoj+i) (essentially, the velocity of 
the system) does not change appreciably over the time 
scale of t c . (This is, basically, the Markov approxima- 
tion.) The result, since the left-hand side of (7) is mul- 
tiplied by ih, is clearly a damping term for (ajcij+i), or, 
by (5), for the on-site velocity Vj = i(a^aj+i — Oj +1 aj): 

d,Vn 

= "27i«j + • • • (10) 

with 

ij = ^((v j+1 vtf) (ii) 
i 



The question now is whether it is possible to extract, 
from the interaction term in (6) something that looks 
like an "external" random potential, as in Eq. (7). That 
this is in fact possible follows if one replaces the bosonic 
field operators hj by hj = Zj + Sj , where Zj is a c- number 
equal to (hj) (the local mean- field), and 5j a zero- average 
operator. Substituting in the interaction term in (6), we 
get 



a](h j+1 - hj)a J+1 = (\z j+1 \ 2 - \zj\ 2 )h]a j+1 

+ a] (z] +1 5 j+ i + z j+1 S j j+1 + 5] +1 5j +1 - z*l 3 - ZjS] - 5}Sj^ h j+1 (12) 

The first term on the right-hand side of (|12|l is a "deterministic" term which can be combined with the first term on 
the right-hand side of Q; indeed, it is the combination of these two terms that yields the Thomas- Fermi profile in 
the ground state when the kinetic energy term (the J term) in Q is negligible. The second term in (|12fl . on the other 
hand, is where we expect the main "noise" to arise. To determine its contribution to the equation of motion for the 
expectation value (ajfij+i), we express the remaining hj operators in terms of Sj, assume expectation values of the 

form ((S' ! ) p 8 q ) vanish unless p — q, and factor terms such as (S^SjSjSj+i) in a standard way, as 2(<5j<5j) (<5j<5_j+i) . The 
result is 

(a](zj +1 S j+1 + z j+ i5]+x + ^j+A+i - z *jh - z i 5 } ~ <%)%+i ) 

= 2 ((s} +1 S j+1 ) (8$)) (z*z j+1 + (6]S j+1 )) + (k,+i| 2 - M 2 ) (Sfa+i) 
= 2(n j+1 - n 3 )(a]a ]+l ) + (\z J+1 \ 2 - \ Zj \ 2 ) (13) 

I 



where the noncondensate density hj = (hj) ~ \zj\ 2 has 
been introduced. The second term on the right-hand side 
of H13|) appears to be a small noise-induced contribution 
to the deterministic part of 112|l . The first term has the 
desired form. We can then replace the expectation value 
of the interaction term in (6) by a deterministic term, 
with which we shall not be concerned any more, and a 
noise-like term 

2U(h j+1 -n j )(a]a j+1 ) (14) 

If the evolution of the hj is sufficiently chaotic, one could 
imagine integrating the Heiscnberg equations many times 
for very slightly different initial conditions and obtaining 
each time a different realization of the "random process" 
hj. Then, if the Markovian condition holds, one can 
follow the same steps as for the external potential Vj in 
Eqs. (7)— (11) above and, by identifying Vj with 2Uhj, 



conclude that, on average, a damping 

2U 2 2U 2 
li = -T2-T c ((hj +1 ~ h f) = -j^r^f 2 ) (15) 

will be observed in this system, for the on-site velocity 
Vj. (For conciseness, we have introduced the notation 
fj = hj + \ — hj). The overall damping of the center of 
mass motion could be estimated by taking a weighted 
average of the jj (although, if the 7^ are very different 
from site to site, the assumption of a single damping 
constant for the center of mass motion may not be a 
very good approximation). 

Our model is, therefore, that the condensate atoms 
are slowed down as they attempt to move through a ran- 
domly fluctuating effective potential created (through the 
interaction term) by the noncondensate atoms, as the lat- 
ter are "shaken" out of equilibrium by the displacement 
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FIG. 2: (a) fi ]+1 -Uj for V = 1E R and j = 0; (b) Fit (dashed 
line) to ln^ 1 [\T{fj (t)] \ 2 ] |, for the first 30 time steps, for 
the fj shown in (a). The slope of this line is taken to be l/r c , 
for this particular value of j. Time is in units of Ti/Er in both 



of the trap. It may be worthwhile, at this point, to go 
over and attempt to justify the various assumptions that 
have been made. 

The interpretation of fij+i — hj as an essentially ran- 
dom variable appears justified from time-dependent HFB 
calculations (about which much more will be said in the 
next section) such as the one illustrated in Fig. 2 for a 
lattice of depth V = 1Er\ the top part shows the time 
trace of hi — ho, and the bottom figure the logarithm 
of the absolute value of its (time-)autocorrelation func- 
tion, with a linear fit showing an approximately exponen- 
tially decaying envelope. It is, however, not quite as clear 
whether the Markov approximation is valid: after all, hj 
is not an external field, but one of the system's dynam- 
ical variables, and it certainly must develop correlations 
and become entangled with other dynamical variables as 
the system evolves. Still, we take this as the simplest 
approximation, and note that, as will be seen in the next 
section, over the range considered, the time scale r c for 



the decay of correlations in n, + i 



(which is, essen- 



tially, h/J) is indeed well-separated from the time scale 
of the damping of the center of mass oscillations. 

Besides the above approximations, we have neglected 
"anomalous averages" such as, e.g., (SjSj), and we have 
used a standard "bosonic" ansatz to factor expectation 
values of products of four operators into expectation val- 
ues of products of two operators. We do this in the spirit 
of all mean-field theories; namely, as something to try 
and see how it works. We certainly do not expect it to 
be a good approximation once (extended) fcrmionization 
becomes important. 

Note that if, instead of using the bosonic asatz, we 
had taken the simplest approach of factoring (dj(hj + \ — 

hj)a,j+i) as ~ (aj(ij + i)(hj + i — hj) (and then separated 
out the condensate part from (rij+i — hj)), the resulting 
"noise" term would have differed from (14) by a factor 
of 2, and hence the damping formula (|15|) would have 
been four times smaller. This may be a reasonable es- 
timate of the possible error involved in our factorization 



assumptions. 

Equation i|15f) does not, by itself, tell us what the ac- 
tual damping is; for that, one needs to know the pa- 
rameters characterizing the strength of the noise, iff), 
and its characteristic correlation time r c . A possible way 
to obtain a very rough order-of-magnitude estimate for 
these quantities has been sketched in |25j |. by rewriting 
Eq. I|15J) in terms of the discrete Fourier transform (mo- 
mentum components) of the hj (the order of magnitude 
of which can be estimated from Fig. 1), and assuming 
that t c should be of the order of magnitude of h/J, since 
this is the "hopping rate," and that is the time scale over 
which one would expect local density fluctuations to de- 
cay. This simple approach does indeed yield the order of 
magnitude of the experimentally-observed damping. 

As we shall show below, using values for (fj) and t c 
derived from HFB calculations in the formula i|15|) does 
lead to very good agreement with the experimentally- 
observed damping. 



IV. HFB CALCULATIONS 

In this section we report on the results of calcula- 
tions using the time dependent Hartree-Fock-Bogoliubov 
(HFB) approximation [22,1221]. The starting point of this 
approximation is the Heiscnberg equation of motion for 
the field operator: 



in — a,- 
dt 3 



(K + ilf + Ua]aj)a 



(16) 



with K the tight binding kinetic energy operator: KAj = 
—J(Aj + i + Aj-i) (here Aj can be any function or oper- 
ator defined at the point j), and flj 2 the external con- 
fining potential, which is quadratic in our system. By 
expressing the field operator, as before, as hj = Zj + 5j, 
replacing this ansatz in Eq. 1|16|) . and treating the cubic 
term in a self-consistent mean field approximation, cou- 
pled equations of motion for the condensate, Zj, and the 
fluctuating field 5j can be obtained: 

ih^-zj = (K + Qj 2 + U(\z 3 \ 2 + 2h j ))z j + Um 3 z*, (17) 

ifl— Sj = (K + Qj 2 + 2Un j )8 j + UmjSl (18) 
dt 



here 



and 



rhj = (SjSj). By using a Bogoliubov transformation 
that expresses the operators Sj in terms of quasiparti- 
cle creation and annihilation operators ctk , &\ and ampli- 

tudes {uj(t)}, {«*(*)}, as Sj = T, k ( u j( t ) & k -vf(t)al), 
one obtains equations of motion for the amplitudes {zj}, 
{u k j (<)} and {vj(t)} known as HFB equations. They de- 
scribe the coupled dynamics of condensate and noncon- 
densate atoms and conserve particle number and energy. 
As coupled, nonlinear equations, they have the potential 
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to describe a wide range of dynamics, including deter- 
ministic chaos. 

A difficulty with the time-independent HFB equations 
is that they violate the Hugenholtz-Pines theorem and 
yield an initial ground state with a depletion that is too 
small p3 | , when compared to the exact numerical results 
in Fig. 1. We have therefore used the Popov approxi- 
mation ps[ (which ignores the anomalous terms fhj) to 
calculate the ground state of the undisplaced trap, but 
then, after displacing the trap, we propagate in time us- 
ing the full HFB equations (without the Popov approx- 
imation), because propagating in time with the Popov 
approximation does not conserve particle number or en- 
ergy. Due to this mixing of approximations (in a sense, 
we are starting from the "wrong" initial state), as well 
as to the intrinsic limitations of the HFB approxima- 
tion, our HFB results must be taken with some caution. 
Nonetheless, one may gain at least some qualitative in- 
sights from them, as illustrated, for instance, in Fig. 3, 
which shows how the noncondensate atoms relax rather 
rapidly (in agreement with the expectation of strongly 
inhibited transport for the high- momentum states), and 
with a substantial amount of noise. 

Our HFB calculations do exhibit damped center of 
mass oscillations for all the values of V in the experi- 
ment, and for sufficiently deep lattices (about V > 4Er 
in our calculations) they even exhibit the overdamped 
relaxation seen in the experiments (i.e., the value of the 
damping exceeds the oscillation frequency), although in 
the experiments this transition to overdamped motion 
was seen already for shallower lattices, between V = 2Er 
and V = 3Er. Quantitatively speaking, the HFB results 
do predict, in general, a smaller damping than is seen ex- 
perimentally for any given lattice depth V, and also, even 
with the Popov approximation, a smaller ground-state 
depletion than the Monte-Carlo calculations in Fig. 1. 
This lack of precise quantitative agreement is not terri- 
bly surprising, given the fact that for all of these systems 
the depiction of the condensate is really not very small 
when compared to the mean field density; hence neglect- 
ing higher powers of the Sj operators cannot be very accu- 
rate. The qualitative agreement, however, suggests that 
the HFB approximation does retain all the physical ingre- 
dients needed to predict the kind of damped oscillations 
seen in the experiments in this regime, 

With all of the above in mind, we have attempted to 
use the results of the HFB calculations to estimate the 
quantities (/J) and r c in the damping formula ll5l in the 
following manner. First, we generate time series for hj{t) 
for all j and for a relatively large number of oscillation 
periods, and we simply average all these values to esti- 
mate (fj(t)). We also calculate the Fourier transform 

fj(u>) = J~{fj) of each time series, and then calculate the 
inverse Fourier transform, J-^ 1 , of the power spectrum 
\fj(u)\ 2 ; by the convolution theorem of Fourier trans- 
forms, this should equal the autocorrelation of fj(t). We 
therefore estimate a correlation time by fitting an expo- 
nential to the decay of the absolute value of J-~ 1 (\fj\ 2 ), 
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FIG. 3: Result of the direct numerical integration of the HFB 
equations for, from top to bottom, V = SEr, V = AEr, and 
V = 5Er. In all cases the position of the center of mass of the 
noncondensate atoms is given by the noisier top trace, that of 
the condensate by the lighter lower trace (small crosses), and 
the total is given by the line drawn with the circles. Time is 
in units of Ti/Er and position is in units of the lattice spacing. 



for relatively short times (of the order of 2>Qh/En). (Rep- 
resentative results are shown in Fig. 2(b).) We obtain in 
this way a (generally different) value of t c and (fj(t)) 
for every lattice site j. The final estimate of the overall 
damping T is obtained by taking a weighted average of 
all the 7j , using the equilibrium density as the weighting 
function. This results in the gray dots in Fig. 4, which 
are to be compared to the experimental data shown as 
the black dots in the same figure. 

We are faced with the somewhat paradoxical result 
that, while the HFB calculations generally underesti- 
mate the damping, the formula (|f 5[) . using HFB values, 
agrees quite well with the experiments and even appears 
to overestimate the damping in places (such as around 
V = 2E R ). 
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FIG. 4: Gray dots (connected by dashed line): the value of 
r calculated from Eq. (3), using the time series that result 
from integrating the HFB equations. Black dots: experimen- 
tal data. 



We do not have, in principle, a reason to doubt the rel- 
ative accuracy of the HFB estimate of the fluctuations' 
correlation time t c (which does turn out to be between 
2H/J and 3h/J for the values of V considered). On the 
other hand, the fact, pointed out above, that the HFB 
calculations predict a noncondcnsatc density lower than 
the true one, suggests that the HFB estimate of (ff(t)) 
may be proportionately low as well. If this is the case, it 
would indicate that the formula H15fl generally overesti- 
mates the damping, perhaps because of the assumption of 
totally uncorrelated condensate and noncondensate fluc- 
tuations that goes into its derivation. The agreement 
with the experiment shown in Fig. 4 would then appear 
to be more precise that it is actually supposed to be. 
Nonetheless, generally speaking, the physical picture in- 
voked in the derivation of the damping formula appears 
to be correct, even if oversimplified in some details (e.g., 
validity of the Markov approximation). 



V. RESULTS FROM GUTZWILLER-ANSATZ 
CALCULATIONS 

It is well-known that an alternative to the HFB calcu- 
lations is provided by a mean-field theory based on the 
Gutzwiller ansatz. While in the HFB method the in- 
teraction term in (1) is treated approximately, and the 
kinetic energy term is treated exactly, yielding a the- 
ory best suited for weakly interacting supcrfluids, the 
Gutzwiller ansatz is equivalent to treating the interac- 
tion term in (1) exactly and approximating the kinetic 
energy term as follows 



(at +1 ) + (at +1 - (ot +1 >) (a,j) + {a 3 - (%}) 



W+i)*aj +a}_ 



(19) 



As a result of this, the Gutzwiller ansatz works best for 
very strongly interacting systems, and it is in fact capable 
of describing qualitatively the Mott transition |30j| , which 
the HFB method cannot do. 
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FIG. 5: Results of Gutzwiller calculations for V = 3Er and 
initial displacements d — 6 (solid line) and d = 8 (dashed 
line). (See text for details.) 



In Eq. I|19fl . the mean field (dj) is obtained self- 
consistently by diagonalization of the resultant effec- 
tive Hamiltonian, to which a chemical potential term 
— /j, ^2 Tij is added in order to get the desired average 
number of particles. Once the initial state has been cal- 
culated, the relevant equations of motion are as given, 
for instance, in fl5L f30| . 

What we find from the Gutzwiller approach is that the 
predicted depletion for the ground state is substantially 
lower than the one calculated numerically in Fig. 1, for 
all except the deepest lattices, and accordingly no appre- 
ciable damping is seen, for the experimental parameters, 
until V = &Er or so. For V = 3-Er a displacement of the 
trap potential by 6 lattice sites fails to give any visible 
damping, but a displacement of 8 lattice sites does give 
substantial damping, as shown in Figure 5(a). 

The other graphs in Figure 5, also for V — 3-Er, high- 
light other interesting features of this transition from un- 
damped to damped motion, which are in rough agree- 
ment with our prior expectations. Figure 5(b) shows the 
fractional depletion as a function of time for the two dis- 
placements d — 6 (solid) and d = 8 (dashed). Although 
the initial depletion is the same (the ground state value), 
the time evolution leads to a depletion that, in the case of 
regular motion, is largest at the times when the conden- 
sate is moving faster. In the case d = 8 one can see the 
depletion initially growing as the speed (the slope of the 
corresponding curve in Fig. 5(a)) increases, and eventu- 
ally becoming rather large, after which damped motion 
follows. Note, for reference, that the Monte-Carlo pre- 
diction from Fig. 1 for this case would be a ground state 
depletion of about 0.5. On the other hand, the smallest 
depletion calculated in Fig. 1, for V = 0.25Er, is 0.24, 



8 



which here would appear to be just large enough to result 
in damped motion. 

Figure 5(c) shows the time dependence of no — n\ (the 
subscript "0" refers to the center of the trap). The reg- 
ular case, for d = 6, is a solid line invisible on the scale 
of the figure (< 0.005 in magnitude). Again, the damp- 
ing appears to be strongly correlated with the site-to-site 
density fluctuations. 

As we did for the HFB calculation, we can extract 
values for (/J) and r c from the time series obtained by 
the Gutzwiller approach, and substitute them in Eq. 1)150. 
The results, for V/Er = 4,5, are actually quite close to 
those obtained from the HFB calculation. For smaller V, 
of course, the calculation would not make sense, since the 
time dependence of nj — hj+\ predicted by the Gutzwiller 
ansatz in this region is always regular, rather than noise- 
like. 



VI. DISCUSSION AND CONCLUSIONS 

To recapitulate, then, we believe that the damping 
in this "intermediate" regime can be explained as aris- 
ing from the large depiction due to the tight trans- 
verse confinement, which leads to the population of high- 
momentum states in the non-quadratic part of the lattice 
dispersion curve. The condensate atoms' motion is then 
damped through their interaction with the random field 
created by these noncondensate atoms when their equi- 
librium state is perturbed. The dramatic growth of Y 
with V illustrated in Fig. 4 arises from several causes: 
first, the depiction increases with lattice depth, as shown 
by Fig. 1 and Eq. J2J); second, the interaction U itself in- 
creases, albeit weakly (as V 1 ^ 4 ); third, the tunneling rate 
J decreases, and the correlation time r c ~ h/J in 115|) 



increases accordingly (the "damping medium" becomes 
more "sluggish"). 

The main limitations of the formula (15) have been 
pointed out when it was derived. Since it only accounts 
for the damping induced by the interactions, it vanishes 
in the U — > limit, even though, as we mentioned in 
the Introduction, a nonintcracting bosonic gas exhibits a 
sort of damping in a lattice, associated with the non- har- 
monic nature of the total potential (see, e.g., @). At the 
other limit point, U — > oo, Eq. (15) predicts an infinite 
damping, which is clearly also not correct. The reason is 
that Eq. (15) is based on a self-consistent factorization 
approximation that, strictly speaking, is only valid in the 
weakly interacting limit. 

We note that, in these regimes were Eq. (15) does not 
apply, previous studies 0, 0, 0, El have shown that 
treatments based on single-particle solutions can provide 
a very accurate description of the damping. On the other 
hand, in the complex intermediate regime where it is not 
possible to use the simplicity of the single-particle so- 
lutions, we have shown that Eq. (15) does manage to 
describe the damping. Moreover, it connects in a sim- 
ple way the damping rate to physical parameters and 
therefore allows a clearer understanding of the physics 
responsible for the dissipative dynamics exhibited by ID 
lattice systems in this regime. 
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